%-------------------------------------------------------------------------------

% This file is part of code_saturne, a general-purpose CFD tool.
%
% Copyright (C) 1998-2025 EDF S.A.
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA 02110-1301, USA.

%-------------------------------------------------------------------------------

\programme{cfxtcl}

\hypertarget{cfxtcl}{}

\vspace{1cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Fonction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Pour le traitement des conditions aux limites, on considère
le système (\ref{Cfbl_Cfxtcl_eq_ref_laminaire_cfxtcl})

\begin{equation}\label{Cfbl_Cfxtcl_eq_ref_laminaire_cfxtcl}
\left\{\begin{array}{l}

\displaystyle\frac{\partial\rho}{\partial t} + \divs(\vect{Q}) = 0 \\
\\
\displaystyle\frac{\partial\vect{Q}}{\partial t}
+ \divv(\vect{u} \otimes \vect{Q}) + \gradv{P}
= \rho \vect{f}_v + \divv(\tens{\Sigma}^v) \\
\\
\displaystyle\frac{\partial E}{\partial t} + \divs( \vect{u} (E+P) )
= \rho\vect{f}_v\cdot\vect{u} + \divs(\tens{\Sigma}^v \vect{u})
- \divs{\,\vect{\Phi}_s} + \rho\Phi_v

\end{array}\right.
\end{equation}

en tant que système hyperbolique portant sur la variable vectorielle
$\vect{W}=\ ^t(\rho,\vect{Q},E)$.

Le système s'écrit alors~:
\begin{equation}\label{Cfbl_Cfxtcl_eq_hyperbolique_cfxtcl}
\displaystyle\frac{\partial \vect{W}}{\partial t}
+ \displaystyle\sum\limits_{i=1}^3
\frac{\partial}{\partial x_i}\vect{F}_i(\vect{W})
= \displaystyle\sum\limits_{i=1}^3
\frac{\partial}{\partial x_i}\vect{F}_i^D(\vect{W},\nabla \vect{W})
+ \vect{\mathcal{S}}
\end{equation}
où les $\vect{F}_i(\vect{W})$ sont les vecteurs flux convectifs
et les $\vect{F}_i^D(\vect{W})$ sont les vecteurs flux diffusifs
dans les trois directions d'espace,
et $\vect{\mathcal{S}}$ est un terme source.

La démarche classique de \CS est adoptée~: on impose les conditions
aux limites en déterminant, pour chaque variable, des valeurs numériques
de bord. Ces valeurs sont calculées de telle façon que, lorsqu'on
les utilise dans les formules standard donnant les flux discrets, on obtienne
les contributions souhaitées au bord.

Pour rendre compte des flux convectifs (aux entrées et aux sorties en particulier),
on fait abstraction des flux diffusifs et des termes
sources pour résoudre un problème de Riemann qui
fournit un vecteur d'état au bord. Celui-ci permet de calculer un flux,
soit directement (par les formules discrètes standard),
soit en appliquant un schéma de Rusanov (schéma de flux décentré).

En paroi, on résout  également, dans certains cas, un problème de
Riemann pour déterminer une pression au bord.

See the \doxygenfile{cfxtcl_8f90.html}{programmers reference of the dedicated subroutine} for further details.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Discrétisation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%=================================
\subsection*{Introduction}
%=================================

%---------------------------------
\subsubsection*{Objectif}
%---------------------------------

On résume ici les différentes conditions aux limites utilisées pour l'algorithme
compressible afin de fournir une vue d'ensemble. Pour atteindre cet objectif,
il est nécessaire de faire référence à  des éléments relatifs
à la discrétisation et au mode d'implantation des conditions aux limites.

Lors de l'implantation, on a cherché à préserver la cohérence avec l'approche
utilisée dans le cadre standard de l'algorithme incompressible de \CS.
Il est donc conseillé d'avoir pris connaissance du mode de traitement des conditions
aux limites incompressibles avant d'aborder les détails de l'algorithme compressible.

Comme pour l'algorithme incompressible, les conditions aux limites sont imposées
par le biais d'une valeur de bord associée à chaque variable. De plus,
pour certaines
frontières (parois à température imposée ou à flux thermique imposé),
on dispose de deux valeurs de bord pour la même variable, l'une d'elles étant dédiée au calcul du
flux diffusif.
Enfin, sur certains types d'entrée et de sortie, on définit également
une valeur du flux convectif au bord.

Comme pour l'algorithme incompressible, l'utilisateur peut définir,
pour chaque face
de bord, des conditions aux limites pour chaque variable, mais on conseille cependant
d'utiliser uniquement les types prédéfinis
décrits ci-après (entrée,
sortie, paroi, symétrie) qui ont l'avantage d'assurer la cohérence entre
les différentes variables et les différentes étapes de calcul.


%---------------------------------
\subsubsection*{Parois}
%---------------------------------

{\bf Pression} : on doit disposer d'une condition pour le calcul du gradient
qui intervient dans l'étape de quantité de mouvement.
On dispose de deux types de condition, au choix de l'utilisateur~:
\begin{itemize}
\item par défaut, la pression imposée au bord est proportionnelle
à la valeur interne (la pression au bord est obtenue comme solution
d'un problème de Riemann sur les équations d'Euler
avec un état miroir~; on distingue les cas de choc et de
détente et, dans le cas d'une détente trop forte, une condition de
Dirichlet homogène est utilisée pour éviter de voir apparaître une
pression négative),
\item si l'utilisateur le souhaite (\var{ICFGRP=1}), le gradient de
pression est imposé à partir du profil de pression hydrostatique.
\end{itemize}
\bigskip

{\bf Vitesse et turbulence}~: traitement standard (voir la documentation des
sous-programmes \fort{cs\_boundary\_conditions}
et \fort{cs\_boundary\_conditions\_set\_coeffs\_turb}).

{\bf Scalaires passifs}~: traitement standard
(flux nul par défaut imposé dans \fort{typecl}).

{\bf Masse volumique}~: traitement standard des scalaires
(flux nul par défaut imposé dans \fort{typecl}).

{\bf Énergie et température\footnote{Le gradient de température est
{\it a priori} inutile, mais peut être requis par l'utilisateur.}}~:
traitement standard des scalaires
(flux nul par défaut imposé dans \fort{typecl}), hormis pour le
calcul du flux diffusif dans le cas de parois à température imposée
ou à flux thermique imposé.

{\bf Flux diffusif pour l'énergie en paroi}~:
l'utilisateur peut choisir (dans \fort{uscfcl}) entre une température de paroi imposée
et un flux thermique diffusif (ou "conductif") imposé.
S'il ne précise rien, on considère que la paroi est adiabatique
(flux thermique diffusif imposé et de valeur nulle).
Dans tous les cas, il faut donc disposer d'un moyen d'imposer le flux diffusif
souhaité. Pour cela, on détermine une valeur de bord pour l'énergie
qui, introduite dans la formule donnant le flux discret, permettra
d'obtenir la contribution attendue
(voir le paragraphe~\ref{Cfbl_Cfxtcl_section_cl_flux_diffusif_energie_cfener}).
Conformément à l'approche classique de \CS, cette valeur est
stockée sous la forme d'un couple de coefficients
(de type \var{COEFAF}, \var{COEFBF}).
Il est important de souligner que cette valeur de bord
ne doit être utilisée que pour le calcul
du flux diffusif~: dans les autres situations pour lesquelles
une valeur de bord de l'énergie ou de la température est requise
(calcul de gradient par exemple), on utilise une condition de flux nul
(traitement standard des scalaires). Pour cela, on dispose d'une
seconde valeur de bord qui est stockée au moyen d'un
couple de coefficients (\var{COEFA}, \var{COEFB}) distinct du précédent.

{\bf Flux convectifs}~: le flux de masse dans la direction normale à la paroi est
pris nul. De ce fait, les flux convectifs seront nuls quelle que soit les valeurs
de bord imposées pour les différentes variables transportées.

%---------------------------------
\subsubsection*{Symétrie}
%---------------------------------

Les conditions appliquées sont les conditions classiques de l'algorithme
incompressible (vitesse normale nulle, flux nul pour les autres variables).

Elles sont imposées dans le sous-programme \fort{typecl} essentiellement. Pour
la pression, la condition de flux nul est imposée dans \fort{cfxtcl}
(au début des développements, on appliquait le même traitement qu'en paroi,
mais une condition de flux nul a été préférée afin de s'affranchir des
problèmes potentiels dans les configurations 2D).


%---------------------------------
\subsubsection*{Entrées et sorties}
%---------------------------------

On obtient, par résolution d'un problème de Riemann au bord, complété par
des relations de thermodynamique (\fort{uscfth}), des valeurs de bord pour toutes
les variables (on suppose qu'en entrée, toutes les composantes de la vitesse
sont fournies~; elles sont supposées nulles par défaut, hormis pour les
entrées à $(\rho,\vect{u})$ imposés, \var{IERUCF}, pour lesquelles il faut
fournir la vitesse explicitement).

Ces valeurs de bord sont utilisées de deux façons~:
\begin{itemize}
\item elles sont utilisées pour calculer les flux convectifs, en faisant
appel au schéma de Rusanov (sauf en sortie supersonique)~; ces flux sont
directement intégrés au second membre des équations à résoudre.
\item elles servent de valeur de Dirichlet dans toutes les autres configurations
pour lesquelles une valeur de bord est requise (calcul de flux diffusif,
calcul de gradient...)
\end{itemize}

Deux cas particuliers~:
\begin{itemize}
\item aux entrées ou sorties pour lesquelles toutes les variables sont imposées
(\var{IESICF}), on utilise une condition de Neumann homogène pour la pression
(hormis pour le calcul du gradient intervenant dans l'équation de la quantité
de mouvement, qui est pris en compte par le flux convectif déterminé par le
schéma de Rusanov). Ce choix est arbitraire (on n'a pas testé le comportement
de l'algorithme si l'on conserve une condition de Dirichlet sur la pression), mais
a été fait en supposant qu'une condition de Neumann homogène serait {\it a
priori} moins déstabilisante, dans la mesure où, pour ce type de frontière,
l'utilisateur peut imposer une valeur de pression très différente de
celle régnant à l'intérieur du domaine (la valeur imposée est utilisée
pour le flux convectif).
\item pour les grandeurs turbulentes et les scalaires utilisateur, si
le flux de masse est entrant et que l'on a fourni
une valeur de Dirichlet (\var{RCODCL(*,*,1)} dans \fort{uscfcl}),
on l'utilise, pour le calcul du flux convectif et du flux diffusif~;
sinon, on utilise une condition de Neumann homogène (le concept
de sortie de type 9 ou 10 est couvert par cette approche).
\end{itemize}



%=================================
\subsection*{Problème de Riemann au bord}
\label{Cfbl_Cfxtcl_section_pb_riemann_cfener}
%=================================

%---------------------------------
\subsubsection*{Introduction}
%---------------------------------

On cherche à obtenir un état au bord,
pour les entrées, les sorties et les parois.

Pour cela, on fait abstraction des flux diffusifs et des sources.
Le système résultant est alors appelé système
d'équations d'Euler. On se place de plus dans un
repère orienté suivant la normale au bord considéré
$(\vect{\tau}_1,\vect{\tau}_2,\vect{n})$
et l'on ne considère que les variations suivant cette normale.
Le système devient donc~:
\begin{equation}\label{Cfbl_Cfxtcl_eq_euler_cfxtcl}
\begin{array}{lllll}
\displaystyle\frac{\partial \vect{W}}{\partial t}
+ \frac{\partial}{\partial n}\vect{F}_n(\vect{W})
= 0
&\text{avec}
& \vect{F}_n(\vect{W})
 = \displaystyle\sum\limits_{i=1}^3 n_i \vect{F}_i(\vect{W})
& \text{et}
& \displaystyle\frac{\partial}{\partial n}
= \displaystyle\sum\limits_{i=1}^3 n_i \frac{\partial}{\partial x_i}
\end{array}
\end{equation}

Pour déterminer les valeurs des variables au bord, on recherche
l'évolution du problème instationnaire suivant,
appelé problème de Riemann~:

\unitlength=1cm
\begin{picture}(20,2.6)
\put(1.5,0){\framebox(12,2.5){}}
\put(7.5,0){\line(0,1){2.5}}
\put(7.5,2.2){bord}
\put(7.5,1.7){\vector(1,0){0.7}}
\put(8.25,1.65){$\vect{n}$}
\multiput(7.5,0)(0,0.5){4}{\line(2,3){.4}}
\put(2,1.2){$\begin{array}{c}
\text{intérieur}\\
\vect{W}_{\,i}\\
\text{état constant dans la cellule $i$}
\end{array}$}
\put(9.5,1.2){$\begin{array}{c}
\text{extérieur}\\
\vect{W}_{\,\infty}\\
\text{état constant}
\end{array}$}
\end{picture}

avec $\vect{W}_{\,\infty}$ dépendant du type de bord et différent
de $\vect{W}_{\,i}$ {\it a priori}.

\vspace{0.3cm}

Pour résoudre ce problème de Riemann, on utilisera les variables
non-conservatives $\widetilde{\vect{W}}=\ ^t(\rho, \vect{u}, P)$
et l'on retrouvera l'énergie grâce à l'équation d'état.

Pour alléger l'écriture, dans le présent
paragraphe~\ref{Cfbl_Cfxtcl_section_pb_riemann_cfener},
on notera aussi $\vect{W}$ le vecteur
$^t(\rho, \vect{u}, P)$ et $\vect{u} = \vect{u}_\tau + u\,\vect{n}$
(en posant $u=\vect{u} \cdot \vect{n}$
et $\vect{u}_\tau = \vect{u} - (\vect{u} \cdot \vect{n})\vect{n}$).

La solution est une suite d'états constants, dont les valeurs
dépendent de $\vect{W}_{\,i}$ et $\vect{W}_{\,\infty}$,
séparés par des ondes se déplaçant à des vitesses données
par les valeurs propres du système $(\lambda_i)_{i=1\ldots 5}$.
On représente les caractéristiques du système sur le schéma suivant~:

\unitlength=1cm
\begin{picture}(20,3.5)
\put(3,0){\vector(1,0){8}}
\put(7,0){\vector(0,1){3}}
\put(11,0.1){$x$}
\put(6.8,2.9){$t$}
\put(7,3.1){bord}
\put(7,2.6){\vector(1,0){0.5}}
\put(7.55,2.65){$\vect{n}$}
\multiput(7,0)(0,0.5){5}{\line(2,3){.4}}
\put(7,0){\line(-1,2){1.4}}
\put(4,2.9){$\lambda_1=u-c$}
\put(7,0){\qbezier[15](0,0)(0.7,1.4)(1.4,2.8)}
\put(8.3,2.9){$\lambda_{2,3,4}=u$}
\put(7,0){\line(2,1){3}}
\put(9.5,1.6){$\lambda_5=u+c$}
\put(5,1.5){$\vect{W}_{\,i}$}
\put(6.3,2){$\vect{W}_{\,1}$}
\put(8.4,1.6){$\vect{W}_{\,2}$}
\put(9.5,0.5){$\vect{W}_{\,\infty}$}
\end{picture}

Comme valeurs des variables au bord, on prendra les valeurs correspondant à
l'état constant qui contient le bord ($\vect{W}_1$ dans l'exemple
précédent).

Il faut remarquer que la solution du problème de Riemann dépend de la
thermodynamique et devra donc être calculée et codée
par l'utilisateur si la thermodynamique n'a pas été prévue
(en version 1.2, la seule
thermodynamique prévue est celle des gaz parfaits).



%---------------------------------
\subsubsection*{En paroi, pour la condition de pression (sans effet de gravité)}
%---------------------------------

Pour les faces de paroi, on définit à l'extérieur du domaine
un état miroir $\vect{W}_{\,\infty}$ par~:
\begin{equation}\label{Cfbl_Cfxtcl_eq_paroi_cfxtcl}
\begin{array}{lllll}
\vect{W}_{\,i} &=&
\left(\begin{array}{l}
\rho_i\\ {\vect{u}_\tau}_i\\ u_i\\ P_i
\end{array}\right)
\qquad
\vect{W}_{\,\infty} &=&
\left(\begin{array}{lll}
\rho_\infty &=& \rho_i\\
{\vect{u}_\tau}_\infty &=& {\vect{u}_\tau}_i\\
u_\infty &=& -u_i\\
P_\infty &=& P_i
\end{array}\right)
\end{array}
\end{equation}

\vspace{0.5cm}

Les solutions dépendent de l'orientation de la vitesse dans la cellule
de bord~:

\vspace{0.5cm}

\begin{enumerate}

\item Si $u_i \leqslant 0$,
la solution est une double détente symétrique.

\unitlength=1cm
\begin{picture}(20,4)
\put(0,0.5){\vector(1,0){8}}
\put(4,0.5){\vector(0,1){3}}
\put(8,0.6){$x$}
\put(3.8,3.4){$t$}
\put(4,3.7){paroi}
\put(4,3.2){\vector(1,0){0.5}}
\put(4.55,3.15){$\vect{n}$}
\multiput(4,0.5)(0,0.5){5}{\line(2,3){.4}}
\put(4,0.5){\line(-1,1){2.5}}
\put(4,0.5){\line(-5,4){2.8}}
\put(4,0.5){\line(-4,5){2}}
\put(0,3.1){$\lambda_1=u-c\ (<0)$}
\put(4.2,2.7){$\lambda_{2,3,4}=0$}
\put(4,0.5){\line(1,1){2.5}}
\put(4,0.5){\line(5,4){2.8}}
\put(4,0.5){\line(4,5){2}}
\put(6.5,3.1){$\lambda_5=u+c\ (>0)$}
\put(1.5,1.5){$\vect{W}_{\,i}$}
\put(3.1,2){$\vect{W}_{\,1}$}
\put(4.5,2){$\vect{W}_{\,2}$}
\put(6,1.5){$\vect{W}_{\,\infty}$}
\put(8.5,2){$\vect{W}_{\,paroi} = \vect{W}_{\,1} = \vect{W}_{\,2}$}
\put(12,2)
{$\left\{\begin{array}{l}
\rho_p = \rho_1 = \rho_2\\
u_p = u_1 = u_2\\
P_p = P_1 = P_2
\end{array}\right.$}
\end{picture}

La conservation de la vitesse tangentielle à travers la 1-onde donne
${\vect{u}_\tau}_p = {\vect{u}_\tau}_i$.
Par des considérations de symétrie on trouve $u_p = 0$.
Puis on obtient $\rho_p$ et $P_p$ en écrivant la conservation
des invariants de Riemann à travers la 1-détente~:
\begin{equation}\label{Cfbl_Cfxtcl_eq_invariants_detente_cfxtcl}
\begin{array}{lll}
\left\{\begin{array}{l}
u_1 + \displaystyle\int_0^{\rho_1} \frac{c}{\rho} d\rho
= u_i + \displaystyle\int_0^{\rho_i} \frac{c}{\rho} d\rho\\
\\
s_1 = s_i
\end{array}\right.
&
\Rightarrow
\left\{\begin{array}{ll}
\displaystyle\int_{\rho_i}^{\rho_1} \frac{c}{\rho} d\rho = u_i
& \Rightarrow \rho_p=\rho_1\\
\\
s(P_1,\rho_1) = s(P_i,\rho_i)
& \Rightarrow P_p=P_1
\end{array}\right.
\end{array}
\end{equation}

\vspace{0.5cm}

\item  Si $u_i > 0$,
la solution est un double choc symétrique.

\unitlength=1cm
\begin{picture}(20,4)
\put(0,0.5){\vector(1,0){8}}
\put(4,0.5){\vector(0,1){3}}
\put(8,0.6){$x$}
\put(3.8,3.4){$t$}
\put(4,3.7){paroi}
\put(4,3.2){\vector(1,0){0.5}}
\put(4.55,3.15){$\vect{n}$}
\multiput(4,0.5)(0,0.5){5}{\line(2,3){.4}}
\put(4,0.5){\line(-1,1){2.5}}
\put(0,3.1){$\lambda_1=u-c\ (<0)$}
\put(4.4,2.7){$\lambda_{2,3,4}=0$}
\put(4,0.5){\line(1,1){2.5}}
\put(6.5,3.1){$\lambda_5=u+c\ (>0)$}
\put(1.5,1.5){$\vect{W}_{\,i}$}
\put(3,2){$\vect{W}_{\,1}$}
\put(4.7,2){$\vect{W}_{\,2}$}
\put(6,1.5){$\vect{W}_{\,\infty}$}
\put(8.5,2){$\vect{W}_{\,paroi} = \vect{W}_{\,1} = \vect{W}_{\,2}$}
\put(12,2)
{$\left\{\begin{array}{l}
\rho_p = \rho_1 = \rho_2\\
u_p = u_1 = u_2\\
P_p = P_1 = P_2
\end{array}\right.$}
\end{picture}

De même que précédemment,
on trouve ${\vect{u}_\tau}_p = {\vect{u}_\tau}_i$ et $u_p = 0$,
puis $\rho_p$ et $P_p$ en écrivant les relations de saut
à travers le 1-choc~:
\begin{equation}\label{Cfbl_Cfxtcl_eq_saut_choc_cfxtcl}
\begin{array}{lll}
\left\{\begin{array}{l}
\rho_1 \rho_i (u_1 - u_i)^2
= (P_1 - P_i)(\rho_1 - \rho_i)\\
\\
2\rho_1 \rho_i (\varepsilon_1 - \varepsilon_i)
= (P_1 + P_i)(\rho_1 - \rho_i)
\end{array}\right.
& \text{avec } \varepsilon = \varepsilon(P,\rho)
&
\Rightarrow
\left\{\begin{array}{l}
\rho_p=\rho_1\\
\\
P_p=P_1
\end{array}\right.
\end{array}
\end{equation}

\bigskip
Pour les gaz parfaits,
avec $M_i = \displaystyle\frac{\vect{u}_i \cdot \vect{n}}{c_i}$
(Nombre de Mach de paroi), on a~:

\begin{itemize}

\item Cas détente ($M_i \leqslant 0$)~:\\
$$
\begin{array}{l}
\left\{\begin{array}{lll}
P_p=0 & \text{si} & 1 + \displaystyle\frac{\gamma-1}{2}M_i<0\\
P_p = P_i \left(1 + \displaystyle\frac{\gamma-1}{2}M_i\right)
^{\frac{2\gamma}{\gamma-1}} & \text{sinon}\\
\end{array}\right.\\
\\
\rho_p=\rho_i \left(\displaystyle\frac{P_p}{P_i}\right)^{\frac{1}{\gamma}}\\
\end{array}
$$

\item Cas choc ($M_i > 0$)~:\\
$$
\begin{array}{l}
P_p = P_i \left(1 + \displaystyle\frac{\gamma(\gamma+1)}{4}M_i^2
+\gamma M_i \displaystyle\sqrt{1+\displaystyle\frac{(\gamma+1)^2}{16}M_i^2}\right)
\\
\\
\rho_p=\rho_i \left(\displaystyle\frac{P_p-P_i}
{P_p-P_i-\rho_i (\vect{u}_i \cdot \vect{n})^2}\right)\\
\end{array}
$$

\end{itemize}

\end{enumerate}

En pratique, le flux convectif normal à la paroi est nul et seule
la condition de pression déterminée ci-dessus est effectivement
utilisée (pour le calcul du gradient sans effet de gravité).

%---------------------------------
\subsubsection*{En sortie}
%---------------------------------

Il existe deux cas de traitement des conditions en sortie,
selon le nombre de Mach normal à la face de bord
($c_i$ est la vitesse du son dans la cellule de bord)~:
$$M_i = \displaystyle\frac{u_i}{c_i}
= \displaystyle\frac{\vect{u}_i \cdot \vect{n}}{c_i}$$

\paragraph{Sortie supersonique (condition ISSPCF de
\fort{uscfcl})~:}
$M_i \geqslant 1 \Rightarrow u_i - c_i \geqslant 0$
\nopagebreak
\linebreak
\unitlength=1cm
\begin{picture}(20,4.5)
\put(0,0.5){\vector(1,0){8}}
\put(4,0.5){\vector(0,1){3}}
\put(8,0.6){$x$}
\put(3.8,3.4){$t$}
\put(4,3.7){bord}
\put(4,3.2){\vector(1,0){0.5}}
\put(4.55,3.15){$\vect{n}$}
\multiput(4,0.5)(0,0.5){5}{\line(2,3){.4}}
\put(4,0.5){\line(1,2){1.4}}
\put(5.3,3.4){$\lambda_1=u-c\ (\geqslant 0)$}
\put(4,0.5){\qbezier[20](0,0)(1.2,1.2)(2.4,2.4)}
\put(6.5,2.9){$\lambda_{2,3,4}=u$}
\put(4,0.5){\line(2,1){3}}
\put(7.1,2){$\lambda_5=u+c$}
\put(2,2){$\vect{W}_{\,i}$}
\put(5.2,2.5){$\vect{W}_{\,1}$}
\put(6,2){$\vect{W}_{\,2}$}
\put(6.5,1){$\vect{W}_{\,\infty}$}
\put(10,2){$\vect{W}_{\,bord} = \vect{W}_{\,i}$}
\end{picture}

Toutes les caractéristiques sont sortantes,
on connaît donc toutes les conditions au bord~:

\begin{equation}
\left\{\begin{array}{l}
\rho_b = \rho_i\\
{\vect{u}_\tau}_b = {\vect{u}_\tau}_i\\
u_b = u_i\\
P_b = P_i
\end{array}\right.
\end{equation}

\paragraph{Sortie subsonique (condition ISOPCF de
\fort{uscfcl})~:}
$0 \leqslant M_i < 1 \Rightarrow (u_i \geqslant 0 \text{ et } u_i - c_i < 0)$

\unitlength=1cm
\begin{picture}(20,4.5)
\put(0,0.5){\vector(1,0){8}}
\put(4,0.5){\vector(0,1){3}}
\put(8,0.6){$x$}
\put(3.8,3.4){$t$}
\put(4,3.7){bord}
\put(4,3.2){\vector(1,0){0.5}}
\put(4.55,3.15){$\vect{n}$}
\multiput(4,0.5)(0,0.5){5}{\line(2,3){.4}}
\put(4,0.5){\line(-1,2){1.4}}
\put(1,3.4){$\lambda_1=u-c\ (<0)$}
\put(4,0.5){\qbezier[15](0,0)(0.7,1.4)(1.4,2.8)}
\put(5.3,3.4){$\lambda_{2,3,4}=u\ (\geqslant 0)$}
\put(4,0.5){\line(2,1){3}}
\put(6.5,2.1){$\lambda_5=u+c$}
\put(2,2){$\vect{W}_{\,i}$}
\put(3.3,2.5){$\vect{W}_{\,1}$}
\put(5.4,2.1){$\vect{W}_{\,2}$}
\put(6.5,1){$\vect{W}_{\,\infty}$}
\put(10,2){$\vect{W}_{\,bord} = \vect{W}_{\,1}$}
\put(12.5,2)
{$\left\{\begin{array}{l}
\rho_b = \rho_1\\
u_b = u_1\\
P_b = P_1
\end{array}\right.$}
\end{picture}

On a une caractéristique entrante,
on doit donc imposer une seule condition  au bord
(en général la pression de sortie $P_{ext}$).

On connaît alors $P_b = P_{ext}$ et ${\vect{u}_\tau}_b = {\vect{u}_\tau}_i$
(par conservation de la vitesse tangentielle à travers la 1-onde).
Pour trouver les inconnues manquantes ($\rho_b$ et $u_b$)
on doit résoudre le passage de la 1-onde~:

\begin{enumerate}

\item Si $P_{ext} \leqslant P_i$,
on a une 1-détente.

On écrit la conservation
des invariants de Riemann à travers la 1-détente~:
\begin{equation}
\begin{array}{lll}
\left\{\begin{array}{l}
s_1 = s_i\\
\\
u_1 + \displaystyle\int_0^{\rho_1} \frac{c}{\rho} d\rho
= u_i + \displaystyle\int_0^{\rho_i} \frac{c}{\rho} d\rho
\end{array}\right.
&
\Rightarrow
\left\{\begin{array}{ll}
s(P_{ext},\rho_1) = s(P_i,\rho_i)
& \Rightarrow \rho_b=\rho_1\\
\\
u_1 = u_i - \displaystyle\int_{\rho_i}^{\rho_1} \frac{c}{\rho} d\rho
& \Rightarrow u_b = u_1
\end{array}\right.
\end{array}
\end{equation}

\item Si $P_{ext} > P_i$,
on a un 1-choc.

On écrit les relations de saut à travers le 1-choc~:
\begin{equation}
\begin{array}{lll}
\left\{\begin{array}{l}
\rho_1 \rho_i (u_1 - u_i)^2
= (P_{ext} - P_i)(\rho_1 - \rho_i)\\
\\
2\rho_1 \rho_i (\varepsilon(P_{ext},\rho_1) - \varepsilon(P_i,\rho_i))
= (P_{ext} + P_i)(\rho_1 - \rho_i)
\end{array}\right.
&
\Rightarrow
\left\{\begin{array}{l}
\rho_b=\rho_1\\
\\
u_b = u_1
\end{array}\right.
\end{array}
\end{equation}

\bigskip
Pour les gaz parfaits, on a~:
\begin{itemize}

\item Cas détente ($P_{ext} \leqslant P_i$)~:\\
$$
\begin{array}{l}
P_b=P_{ext}\\
\\
\rho_b=\rho_i \left(\displaystyle\frac{P_{ext}}{P_i}\right)^{\frac{1}{\gamma}}
\end{array}
$$

\item Cas choc ($P_{ext} > P_i$)~:\\
$$
\begin{array}{l}
P_b=P_{ext}\\
\\
\rho_b=\rho_i \left(\displaystyle\frac{P_{ext}-P_i}{P_{ext}-P_i-\rho_i
(\vect{u}_i \cdot \vect{n} - \vect{u}_b \cdot \vect{n})^2}\right)
= \rho_i \left(\displaystyle\frac{(\gamma+1)P_{ext}+(\gamma-1)P_i}
{(\gamma-1)P_{ext}+(\gamma+1)P_i}\right)\\
\end{array}
$$

\end{itemize}

\end{enumerate}

La valeur de la masse volumique au bord intervient en particulier
dans le flux de masse.


%---------------------------------
\subsubsection*{En entrée}
%---------------------------------

L'utilisateur impose les valeurs qu'il souhaite pour les variables
en entrée~:
$$
\begin{array}{lllll}
\vect{W}_{\,ext} &=&
\left(\begin{array}{l}
\rho_{ext}\\ {\vect{u}_\tau}_{ext}\\ u_{ext}\\ P_{ext}
\end{array}\right)
\end{array}
$$

De même que précédemment, il existe deux cas de traitement
des conditions en entrée,
pilotés par le nombre de Mach entrant, normalement à la face de bord
(avec $c_{ext}$ la vitesse du son en entrée)~:
$$M_{ext} = \displaystyle\frac{u_{ext}}{c_{ext}}
= \displaystyle\frac{\vect{u}_{ext} \cdot \vect{n}}{c_{ext}}$$

\paragraph{Entrée supersonique (condition IESICF de
\fort{uscfcl})~:}
$M_{ext} \leqslant -1 \Rightarrow u_{ext} + c_{ext} \leqslant 0$

\unitlength=1cm
\begin{picture}(20,4.5)
\put(0,0.5){\vector(1,0){8}}
\put(4,0.5){\vector(0,1){3}}
\put(8,0.6){$x$}
\put(3.8,3.4){$t$}
\put(4,3.7){bord}
\put(4,3.2){\vector(1,0){0.5}}
\put(4.55,3.15){$\vect{n}$}
\multiput(4,0.5)(0,0.5){5}{\line(2,3){.4}}
\put(4,0.5){\line(-2,1){3}}
\put(0,2.1){$\lambda_1=u-c$}
\put(4,0.5){\qbezier[20](0,0)(-1.2,1.2)(-2.4,2.4)}
\put(0,2.9){$\lambda_{2,3,4}=u$}
\put(4,0.5){\line(-1,2){1.4}}
\put(1,3.4){$\lambda_5=u+c\ (\leqslant 0)$}
\put(1.5,1){$\vect{W}_{\,i}$}
\put(2,1.7){$\vect{W}_{\,1}$}
\put(2.2,2.5){$\vect{W}_{\,2}$}
\put(6,2){$\vect{W}_{\,\infty}$}
\put(10,2){$\vect{W}_{\,bord} = \vect{W}_{\,\infty} = \vect{W}_{\,ext}$}
\end{picture}

Toutes les caractéristiques sont entrantes,
toutes les conditions au bord sont donc imposées par l'utilisateur.

\begin{equation}
\left\{\begin{array}{l}
\rho_b = \rho_{ext}\\
{\vect{u}_\tau}_b = {\vect{u}_\tau}_{ext}\\
u_b = u_{ext}\\
P_b = P_{ext}
\end{array}\right.
\end{equation}


\paragraph{Entrée subsonique (condition IERUCF de
\fort{uscfcl})~: }
$$-1 < M_{ext} \leqslant 0
\Rightarrow (u_{ext} \leqslant 0 \text{ et } u_{ext} + c_{ext} > 0)$$


\unitlength=1cm
\begin{picture}(20,4.5)
\put(0,0.5){\vector(1,0){8}}
\put(4,0.5){\vector(0,1){3}}
\put(8,0.6){$x$}
\put(3.8,3.4){$t$}
\put(4,3.7){bord}
\put(4,3.2){\vector(1,0){0.5}}
\put(4.55,3.15){$\vect{n}$}
\multiput(4,0.5)(0,0.5){5}{\line(2,3){.4}}
\put(4,0.5){\line(-2,1){3}}
\put(0,2.1){$\lambda_1=u-c$}
\put(4,0.5){\qbezier[15](0,0)(-0.7,1.4)(-1.4,2.8)}
\put(1.1,3.4){$\lambda_{2,3,4}=u\ (\leqslant 0)$}
\put(4,0.5){\line(1,2){1.4}}
\put(5.3,3.4){$\lambda_5=u+c\ (>0)$}
\put(1.5,1){$\vect{W}_{\,i}$}
\put(2,2.1){$\vect{W}_{\,1}$}
\put(3.3,2.5){$\vect{W}_{\,2}$}
\put(6,2){$\vect{W}_{\,\infty}$}
\put(10,2){$\vect{W}_{\,bord} = \vect{W}_{\,2}$}
\put(12.5,2)
{$\left\{\begin{array}{l}
\rho_b = \rho_2\\
u_b = u_2\\
P_b = P_2
\end{array}\right.$}
\end{picture}


On a une caractéristique sortante.
L'utilisateur doit donc laisser un degré de liberté.

En général, on impose le flux de masse en entrée, donc $\rho_{ext}$
et $u_{ext}$, et l'on calcule la pression au bord en résolvant
le passage des 1$\sim$4-ondes.
On connaît aussi ${\vect{u}_\tau}_b = {\vect{u}_\tau}_{ext}$,
par conservation de la vitesse tangentielle à travers la 5-onde.

\begin{enumerate}

\item Si $u_{ext} \geqslant u_i$,
on a une 1-détente.

On écrit la conservation
des invariants de Riemann à travers la 1-détente
et la conservation de la vitesse et de la pression à travers le contact~:
\begin{equation}
\begin{array}{lll}
\begin{array}{l}
\left\{\begin{array}{l}
u_1 + \displaystyle\int_0^{\rho_1} \frac{c}{\rho} d\rho
= u_i + \displaystyle\int_0^{\rho_i} \frac{c}{\rho} d\rho\\
\\
s_1 = s_i
\end{array}\right.\\
\\
\left\{\begin{array}{l}
u_1 = u_2 = u_{ext}\\
\\
P_1 = P_2
\end{array}\right.
\end{array}
&
\Rightarrow
\left\{\begin{array}{ll}
\displaystyle\int_{\rho_i}^{\rho_1} \frac{c}{\rho} d\rho
= u_i - u_{ext}
& \Rightarrow \rho_1\\
\\
s(P_2,\rho_1) = s(P_i,\rho_i)
& \Rightarrow P_b = P_2
\end{array}\right.
\end{array}
\end{equation}

\item Si $u_{ext} < u_i$,
on a un 1-choc.

On écrit les relations de saut à travers le 1-choc
et la conservation de la vitesse et de la pression à travers le contact~:
\begin{equation}
\begin{array}{lll}
\begin{array}{l}
\left\{\begin{array}{l}
\rho_1 \rho_i (u_1 - u_i)^2
= (P_1 - P_i)(\rho_1 - \rho_i)\\
\\
2\rho_1 \rho_i (\varepsilon_1 - \varepsilon_i)
= (P_1 + P_i)(\rho_1 - \rho_i)\\
\\
\varepsilon = \varepsilon(P,\rho)
\end{array}\right.\\
\\
\left\{\begin{array}{l}
u_1 = u_2 = u_{ext}\\
\\
P_1 = P_2
\end{array}\right.
\end{array}
&
\Rightarrow
\left\{\begin{array}{l}
\rho_1\\
\\
P_b = P_2
\end{array}\right.
\end{array}
\end{equation}

\bigskip
Pour les gaz parfaits, on a~:

\begin{itemize}

\item Cas détente ($\delta M \leqslant 0$)~:\\
$$
\begin{array}{l}
\left\{\begin{array}{lll}
P_b=0 & \text{si} & 1 + \displaystyle\frac{\gamma-1}{2}\delta M<0\\
P_b = P_i \left(1 + \displaystyle\frac{\gamma-1}{2}\delta M\right)
^{\frac{2\gamma}{\gamma-1}} & \text{sinon}\\
\end{array}\right.\\
\\
\rho_b= \rho_{ext}\\
\end{array}
$$

\item Cas choc ($\delta M > 0$)~:\\
$$
\begin{array}{l}
P_b = P_i \left(1 + \displaystyle\frac{\gamma(\gamma+1)}{4}\delta M^2
+\gamma \delta M \displaystyle\sqrt{1+\displaystyle\frac{(\gamma+1)^2}{16}\delta M^2}\right)
\\
\\
\rho_b=\rho_{ext}\\
\end{array}
$$

\end{itemize}

\end{enumerate}


%=================================
\subsection*{Condition de pression en paroi avec effets de gravité}
%=================================

Le problème de Riemann considéré précédemment
ne prend pas en compte les effets de la gravité.
Or, dans certains cas, si l'on ne prend pas en compte le gradient de
pression ``hydrostatique'',  on peut obtenir une solution erronée
(en particulier, par exemple,
on peut créer une source de quantité de mouvement
non physique dans un milieu initialement au repos).

Écrivons l'équilibre local dans la maille de bord~:
\begin{equation}\label{Cfbl_Cfxtcl_eq_equilibre_local_cfxtcl}
\gradv{P} = \rho \vect{g}
\end{equation}

Pour simplifier la résolution, on peut utiliser la formulation
de (\ref{Cfbl_Cfxtcl_eq_equilibre_local_cfxtcl}) en incompressible
(c'est cette approche qui a été adoptée dans \CS)~:
\begin{equation}\label{Cfbl_Cfxtcl_eq_equilibre_incompressible_cfxtcl}
\begin{array}{lll}
\left(\gradv{P}\right)_i = \rho_i \vect{g}
& \text{ce qui donne}
& P_{paroi} = P_i + \rho_i \vect{g} \cdot (\vect{x}_{paroi} - \vect{x}_i)
\end{array}
\end{equation}


Une autre approche (dépendante de l'équation d'état)
consiste à résoudre l'équilibre local avec la formulation
compressible (\ref{Cfbl_Cfxtcl_eq_equilibre_local_cfxtcl}), en supposant de plus que
la maille est isentropique~:
\begin{equation}
\left\{\begin{array}{lll}
\gradv{P} = \rho \vect{g}\\
\\
P = P(\rho,s_i)
\end{array}\right.
\end{equation}
Ce qui donne, pour un gaz parfait~:
\begin{equation}
\label{Cfbl_Cfxtcl_eq_equilibre_compressible_cfxtcl}
 P_{paroi} = P_i \left( 1+ \displaystyle\frac{\gamma -1}{\gamma}
\displaystyle\frac{\rho_i}{P_i} \vect{g} \cdot (\vect{x}_{paroi} - \vect{x}_i)
\right)^{\frac{\gamma}{\gamma -1}}
\end{equation}

\paragraph{Remarque~:}
la formule issue de l'incompressible (\ref{Cfbl_Cfxtcl_eq_equilibre_incompressible_cfxtcl})
est une linéarisation de la formule (\ref{Cfbl_Cfxtcl_eq_equilibre_compressible_cfxtcl}).
Dans les cas courants elle s'éloigne très peu de la formule exacte.
Dans des conditions extrêmes,
si l'on considère par exemple
de l'air à $1000K$ et $10bar$, avec une accélération
de la pesanteur $g=1000m/s^2$ et une différence de hauteur entre
le centre de la cellule et le centre de la face de bord de $10m$,
l'expression (\ref{Cfbl_Cfxtcl_eq_equilibre_compressible_cfxtcl}) donne $P_{paroi} = 1034640,4Pa$
et l'expression (\ref{Cfbl_Cfxtcl_eq_equilibre_incompressible_cfxtcl}) donne $P_{paroi} = 1034644,7Pa$,
soit une différence relative de moins de $0,001\%$.
On voit aussi que la différence entre la pression calculée au centre
de la cellule et celle calculée au bord est de l'ordre de~$3\%$.

%=================================
\subsection*{Schéma de Rusanov pour le calcul de flux convectifs au bord}
%=================================


%---------------------------------
\subsubsection*{Introduction}
%---------------------------------

Le schéma de Rusanov est utilisé pour certains types de conditions aux
limites afin de passer du vecteur d'état calculé au bord comme indiqué
précédemment (solution du problème de Riemann) à un flux convectif de
bord (pour la masse, la quantité de
mouvement et l'énergie). L'utilisation de ce schéma (décentré amont)
permet de gagner en  stabilité.

Le schéma de Rusanov est appliqué aux frontières auxquelles on considère
qu'il est le plus probable de rencontrer des conditions en accord imparfait
avec l'état régnant dans le domaine, conditions qui sont donc susceptibles de
déstabiliser le calcul~: il s'agit des entrées et des sorties (frontières
de type IESICF, ISOPCF, IERUCF, IEQHCF). En sortie
supersonique (ISSPCF) cependant, le schéma de Rusanov est inutile et
n'est donc pas appliqué~:
en effet, pour ce type de frontière, l'état imposé au bord est exactement
l'état amont et le décentrement du schéma de Rusanov n'apporterait donc
rien.

%---------------------------------
\subsubsection*{Principe}
%---------------------------------

Pour le calcul du flux décentré de Rusanov, on considère
le système hyperbolique
constitué des seuls termes convectifs issus
des équations de masse, quantité de mouvement et énergie. Ce
système est écrit, par changement de variable, en non conservatif
(on utilise la relation
$\displaystyle P=\frac{\rho\varepsilon}{\gamma-1}$ et
on note $u_\xi$ les composantes de $\vect{u}$)~:

\begin{equation}
\left\{\begin{array}{lllll}
\displaystyle\frac{\partial\rho}{\partial t}
&+&\rho\divv{\,\vect{u}} + \vect{u}\,\grad{\,\rho}&=& 0 \\
\displaystyle\frac{\partial u_\xi}{\partial t}
&+& \vect{u}\,\grad{u_\xi}+\displaystyle\frac{1}{\rho}\,\frac{\partial
P}{\partial \xi} &=& 0 \\
\displaystyle\frac{\partial P}{\partial t}
&+&\gamma\,P\,\dive{\vect{u}}+\vect{u}\,\grad{P}&=& 0
\end{array}\right.
\end{equation}

En notant le vecteur d'état $\vect{W}= (\rho,\vect{u},P)^t$,
ce système est noté~:
\begin{equation}
\displaystyle\frac{\partial \vect{W}}{\partial t} +\dive{\,\vect{F}(\vect{W})} = 0
\end{equation}

Avec $\delta\,\vect{W}$ l'incrément temporel du vecteur d'état, $\vect{n}$ la
normale à une face, $ij$ la face interne partagée par les cellules $i$ et
$j$ et $ik$ la face de bord $k$ associée à la cellule $i$,
la discrétisation spatiale conduit à~:
\begin{equation}
\displaystyle\frac{|\Omega_i|}{\Delta t}\delta\,\vect{W}_i
+\sum\limits_{j\in\,Vois(i)}\int_{S_{ij}} \vect{F}(\vect{W})\,\vect{n}\,dS
+\sum\limits_{k\in {\gamma_b(i)}}\int_{S_{\,{b}_{ik}}} \vect{F}(\vect{W})\,\vect{n}\,dS
=0
\end{equation}

Sur une face de bord donnée,
on applique le schéma de Rusanov pour calculer le flux
comme suit~:
\begin{equation}
\frac{1}{|S_{\,{b}_{ik}}|}\int_{S_{\,{b}_{ik}}} \vect{F}(\vect{W})\,\vect{n}\,dS
=\frac{1}{2}\left(\vect{F}(\vect{W}_i)+\vect{F}(\vect{W}_{\,{b}_{ik}})\right)\cdot\vect{n}_{\,{b}_{ik}}
-\frac{1}{2}\rho_{rus\,{b}_{ik}}\left(\vect{W}_{\,{b}_{ik}}-\vect{W}_i\right)=\vect{F}_{rus\,{b}_{ik}}(\vect{W})
\end{equation}

Dans cette relation, $\vect{W}_{\,{b}_{ik}}$ est  le vecteur d'état
$\vect{W}_{\infty}$, connu au bord (tel
qu'il résulte de la résolution du problème de Riemann au bord
présentée plus haut pour chaque type de frontière considéré).

%---------------------------------
\subsubsection*{Paramètre de décentrement $\rho_{rus\,{b}_{ik}}$}
%---------------------------------

Pour chaque face de bord, le scalaire $\rho_{rus\,{b}_{ik}}$ est la
plus grande valeur du rayon spectral de la matrice jacobienne
$\displaystyle\frac{\partial\,\vect{F}_n(\vect{W})}{\partial \vect{W}}$
obtenu pour les vecteurs d'état $\vect{W}_i$ et $\vect{W}_{\,{b}_{ik}}$.

$\vect{F}_n$ est la composante du
flux $\vect{F}$ dans la direction de la normale à la face de bord,
$\vect{n}_{\,{b}_{ik}}$. Utiliser $\vect{F}_n$
pour la détermination du
paramètre de décentrement $\rho_{rus\,{b}_{ik}}$
relève d'une approche classique qui consiste
à remplacer le système tridimensionnel
initial par le système unidimensionnel projeté dans la direction
normale à la face, en négligeant les variations du vecteur d'état
$\vect{W}$ dans la direction tangeante à la face~:
\begin{equation}
\displaystyle\frac{\partial \vect{W}}{\partial t} +\frac{\partial\,\vect{F}_n(\vect{W})}{\partial
\vect{W}}\,\frac{\partial \vect{W}}{\partial n} = 0
\end{equation}

De manière plus explicite, si l'on se place dans un repère de calcul ayant
$\vect{n}_{\,{b}_{ik}}$ comme  vecteur de base, et si l'on note $u$ la
composante de vitesse associée, le système est le suivant (les équations
portant sur les composantes transverses de la vitesse sont découplées,
associées à la valeur propre $u$, comme le serait un scalaire simplement
convecté et ne sont pas écrites ci-après)~:
\begin{equation}
\left\{\begin{array}{lllll}
\displaystyle\frac{\partial\rho}{\partial t}
&+&\displaystyle\rho\frac{\partial\,u}{\partial\,n} + u\,\frac{\partial\,\rho}{\partial\,n}&=& 0 \\
\displaystyle\frac{\partial u}{\partial t}
&+&\displaystyle u\,\frac{\partial\,u}{\partial\,n}+\frac{1}{\rho}\,\frac{\partial
P}{\partial n} &=& 0 \\
\displaystyle\frac{\partial P}{\partial t}
&+&\displaystyle\gamma\,P\,\frac{\partial\,u}{\partial\,n}+u\,\frac{\partial\,P}{\partial\,n}&=& 0
\end{array}\right.
\end{equation}

La matrice jacobienne associée est donc~:
\begin{equation}
\left(\begin{array}{lll}
\displaystyle u & \rho                 & 0                                \\
\displaystyle 0 & u                    & \displaystyle\frac{1}{\rho}        \\
\displaystyle 0 & \gamma\, P        & 0                                 \\
\end{array}\right)
\end{equation}

Les valeurs propres sont $u$ et $\displaystyle\,u\pm c$ (avec
$c=\sqrt\frac{\gamma\,P}{\rho}$). Le rayon spectral est donc
$|u|+c$ et le paramètre de décentrement s'en déduit~:
\begin{equation}
\rho_{rus\,{b}_{ik}} = max\left(|u_i|+c_i,|u_{{b}_{ik}}|+c_{{b}_{ik}}\right)
\end{equation}


%---------------------------------
\subsubsection*{Expression des flux convectifs}
%---------------------------------

Les flux convectifs calculés par le schéma de Rusanov
pour les variables masse, quantité de mouvement
et énergie représentent donc la discrétisation des termes suivants~:
\begin{equation}
\left\{\begin{array}{l}
\displaystyle\dive(\vect{Q})\\
\displaystyle\divv(\vect{u}\otimes\vect{Q})+\grad\,P\\
\displaystyle\dive\left(\vect{Q}\,(e+\frac{P}{\rho})\right)
\end{array}\right.
\end{equation}

Pour une face de bord $ik$ adjacente à la cellule $i$ et
avec la valeur précédente de $\rho_{rus\,{b}_{ik}}$, on a~:
\begin{equation}
\left\{\begin{array}{lll}
\displaystyle\int_{S_{\,{b}_{ik}}}\vect{Q}\cdot\vect{n}\,dS
&=&
\displaystyle\frac{1}{2}\left(
(\vect{Q}_i+\vect{Q}_{\,{b}_{ik}})\cdot\vect{n}_{\,{b}_{ik}}\right)\,S_{\,{b}_{ik}}\\
&&\displaystyle \qquad\qquad\qquad\qquad\qquad\qquad
-\frac{1}{2}\,\rho_{rus\,{b}_{ik}}
\left(\rho_{\,{b}_{ik}}-\rho_i \right)\,S_{\,{b}_{ik}}\\
%
\displaystyle\int_{S_{\,{b}_{ik}}}(\vect{u}\otimes\vect{Q}+\grad\,P)\cdot\vect{n}\,dS
&=&
\displaystyle\frac{1}{2}\left(
 \vect{u}_i(\vect{Q}_i\cdot\vect{n}_{\,{b}_{ik}})
+P_i\,\vect{n}_{\,{b}_{ik}}
+\vect{u}_{\,{b}_{ik}}(\vect{Q}_{\,{b}_{ik}}\cdot\vect{n}_{\,{b}_{ik}})
+P_{\,{b}_{ik}}\,\vect{n}_{\,{b}_{ik}}\right)\,S_{\,{b}_{ik}}\\
&&\displaystyle \qquad\qquad\qquad\qquad\qquad\qquad
-\frac{1}{2}\,\rho_{rus\,{b}_{ik}}
\left(\vect{Q}_{\,{b}_{ik}}-\vect{Q}_i \right)\,S_{\,{b}_{ik}}\\
%
\displaystyle\int_{S_{\,{b}_{ik}}}(e+\frac{P}{\rho})\,\vect{Q}\cdot\vect{n}\,dS
&=&
\displaystyle\frac{1}{2}\left(
(e_i+\frac{P_i}{\rho_i})\,(\vect{Q}_i\cdot\vect{n}_{\,{b}_{ik}})
+(e_{\,{b}_{ik}}+\frac{P_{\,{b}_{ik}}}{\rho_{\,{b}_{ik}}})(\vect{Q}_{\,{b}_{ik}}\cdot\vect{n}_{\,{b}_{ik}})
\right)\,S_{\,{b}_{ik}}\\
&&\displaystyle \qquad\qquad\qquad\qquad\qquad\qquad
-\frac{1}{2}\,\rho_{rus\,{b}_{ik}}
\left(\rho_{\,{b}_{ik}}\,e_{\,{b}_{ik}}-\rho_i\,e_i \right)\,S_{\,{b}_{ik}}\\
\end{array}\right.
\end{equation}



%=================================
\subsection*{Conditions aux limites pour le flux diffusif d'énergie}
\label{Cfbl_Cfxtcl_section_cl_flux_diffusif_energie_cfener}
%=================================

%---------------------------------
\subsubsection*{Rappel}
%---------------------------------

Pour le flux de diffusion d'énergie, les conditions aux limites sont
imposées de manière similaire à ce qui est décrit dans
la documentation de \fort{cs\_boundary\_conditions\_set\_coeffs\_turb} et de
\fort{cs\_boundary\_conditions}. La figure~(\ref{Cfbl_Cfxtcl_fig_flux_cfxtcl}) rappelle quelques notations
usuelles et l'équation~(\ref{Cfbl_Cfxtcl_eq_flux_cfxtcl}) traduit la conservation du flux
normal au bord pour la variable $f$.

\begin{figure}[htp]
\centerline{\includegraphics[height=7cm]{fluxbord}}
\caption{\label{Cfbl_Cfxtcl_fig_flux_cfxtcl}Cellule de bord.}
\end{figure}

\begin{equation}\label{Cfbl_Cfxtcl_eq_flux_cfxtcl}
\begin{array}{l}
    \underbrace{h_{int}(f_{b,int}-f_{I'})}_{\phi_{int}}
  = \underbrace{h_{b}(f_{b,ext}-f_{I'})}_{\phi_{b}}
  = \left\{\begin{array}{ll}
    \underbrace{h_{imp,ext}(f_{imp,ext}-f_{b,ext})}_{\phi_{\text{\it réel
imposé}}} &\text{(condition de Dirichlet)}\\
    \underbrace{\phi_{\text{\it imp,ext}}}_{\phi_{\text{\it réel imposé}}}
            &\text{(condition de Neumann)}
           \end{array}\right.
\end{array}
\end{equation}


L'équation~(\ref{Cfbl_Cfxtcl_eq_fbint_cfxtcl}) rappelle la formulation des
conditions aux limites pour une variable $f$.
\begin{equation}\label{Cfbl_Cfxtcl_eq_fbint_cfxtcl}
f_{b,int}
  = \left\{\begin{array}{cccccl}
    \displaystyle\frac{h_{imp,ext}}{h_{int}+h_r h_{imp,ext} }&f_{imp,ext}&+&
    \displaystyle\frac{h_{int}+h_{imp,ext}(h_r-1)}{h_{int}+h_r h_{imp,ext} }&f_{I'}
                         &\text{(condition de Dirichlet)}\\
    \displaystyle\frac{1}{h_{int}}&\phi_{\text{\it imp,ext}}&+&
    \ &f_{I'}
            &\text{(condition de Neumann)}
           \end{array}\right.
\end{equation}

Les coefficients d'échange sont définis comme suit\footnote{On rappelle que, comme
dans \fort{cs\_boundary\_conditions}, $\alpha$ désigne $\lambda+C_p\,\frac{\mu_t}{\sigma_t}$
si $f$ est la température,
$\frac{\lambda}{C_p}+\frac{\mu_t}{\sigma_t}$ si $f$ représente l'enthalpie.
Le coefficient $C$ représente $C_p$ pour la température et vaut
$1$ pour l'enthalpie. La grandeur adimensionnelle $f^+$ est obtenue par
application d'un principe de similitude en paroi~: pour la température,
elle dépend du nombre de
Prandlt moléculaire, du nombre de Prandtl turbulent et de la distance adimensionnelle à la paroi $y^+$ dans la cellule de bord.}~:
\begin{equation}
\left\{\begin{array}{lll}
h_{int}&=&\displaystyle\frac{\alpha}{\overline{I'F}}\\
h_r&=&\displaystyle\frac{h_{int}}{h_{b}} \\
h_b&=&\displaystyle\frac{\phi_b}{f_{b,ext}-f_{I'}}=\frac{\rho\,C\,u_k}{f^+_{I'}}
\end{array}\right.
\end{equation}

Dans \CS, on note les conditions aux limites de manière générale sous
la forme suivante~:
\begin{equation}
f_{b,int}=A_b + B_b\,f_{I'}
\end{equation}
avec $A_b$ et $B_b$ définis selon le type des conditions~:
\begin{equation}
\begin{array}{c}
\text{Dirichlet}\left\{\begin{array}{ll}
    A_b = &\displaystyle\frac{h_{imp,ext}}{h_{int}+h_r h_{imp,ext} } f_{imp,ext}\\
    B_b = &\displaystyle\frac{h_{int}+h_{imp,ext}(h_r-1)}{h_{int}+h_r h_{imp,ext} }
                  \end{array}\right.
\text{\ \  Neumann}\left\{\begin{array}{ll}
    A_b = &\displaystyle\frac{1}{h_{int}}\phi_{\text{\it imp,ext}}\\
    B_b = &1
                  \end{array}\right.
\end{array}
\end{equation}

%---------------------------------
\subsubsection*{Flux diffusif d'énergie}
%---------------------------------

Dans le module compressible, on résout une équation sur l'énergie, qui s'écrit, si
l'on excepte tous les termes hormis le flux de diffusion et le terme
instationnaire, pour faciliter la présentation~:

\begin{equation}
\begin{array}{lll}
\displaystyle\frac{\partial \rho e}{\partial t} &=& - \dive{\,\vect{\Phi}_s}\\
&=& \displaystyle\dive{(K\,\grad{T})} \text{\ \ avec \ \ } K=\lambda+C_p\,\frac{\mu_t}{\sigma_t}\\
&=& \displaystyle\dive{\left(K\,\grad{\frac{e-\frac{1}{2}\,u^2-\varepsilon_{sup}}{C_v}}\right)} \\
&=& \displaystyle\dive{\left(\frac{K}{C_v}\,\grad{(e-\frac{1}{2}\,u^2-\varepsilon_{sup})}\right)} \text{\ \
si \ \ } C_v \text{\ est constant}\\
&=& \displaystyle\dive{\left(\frac{K}{C_v}\,\grad\,e\right)}
-\dive{\left(\frac{K}{C_v}\,\grad{(\frac{1}{2}\,u^2+\varepsilon_{sup})}\right)} \\

\end{array}
\end{equation}

La décomposition en $e$ et $\frac{1}{2}\,u^2+\varepsilon_{sup}$ est purement
mathématique (elle résulte du fait que l'on résout en énergie alors que
le flux thermique s'exprime en fonction de la température). Aussi,  pour imposer un
flux de bord ou une température de bord (ce qui revient au même puisque l'on
impose toujours finalement la conservation du flux normal), on {\it choisit}
de reporter la totalité de la condition à la limite sur le terme
$\displaystyle\frac{K}{C_v}\,\grad\,e$
et donc d'annuler le flux associé au terme
$\displaystyle\frac{K}{C_v}\,\grad{(\frac{1}{2}\,u^2+\varepsilon_{sup})}$
(en pratique, pour l'annuler, on se contente de ne pas l'ajouter
au second membre de l'équation). Conformément à l'approche retenue dans \CS et
rappelée précédemment, on déterminera donc une valeur de bord {\it
fictive} de l'énergie qui permette de reconstruire le flux diffusif total
attendu à partir
de la discrétisation du seul terme $\displaystyle\frac{K}{C_v}\,\grad\,e$.

Remarque : dans la version 1.2.0,
on utilise $\displaystyle
\frac{K}{C_v}=\left(\frac{\lambda}{C_v}+\frac{\mu_t}{\sigma_t}\right)$, à
partir de 1.2.1, on utilise la valeur  $\displaystyle
\frac{K}{C_v}=\left(\frac{\lambda}{C_v}+\frac{C_p}{C_v}\frac{\mu_t}{\sigma_t}\right)$.
On notera que le nombre de Prandtl turbulent $\sigma_t$ est associé à la variable
résolue et peut être fixé par l'utilisateur.


%---------------------------------
\subsubsection*{Condition de Neumann}
%---------------------------------

La conservation du flux s'écrit~:

\begin{equation}
    \underbrace{h_{int}(e_{b,int}-e_{I'})}_{\phi_{int}}
    =\underbrace{\phi_{\text{\it imp,ext}}}_{\phi_{\text{\it réel imposé}}}
\end{equation}

On a donc dans ce cas~:
\begin{equation}
\left\{\begin{array}{lll}
  A_b &= &\displaystyle\frac{1}{h_{int}}\phi_{\text{\it imp,ext}}\\
  B_b &= &1
\end{array}\right.
\end{equation}


%---------------------------------
\subsubsection*{Condition de Dirichlet}
%---------------------------------

On suppose que la condition de Dirichlet porte sur la température $T_{b,ext}$.


La conservation du flux s'écrit~:
\begin{equation}\label{Cfbl_Cfxtcl_eq_conservation_flux_cfxtcl}
    \underbrace{h_{int}(e_{b,int}-e_{I'})}_{\phi_{int}\text{\ (forme numérique
du flux)}}
  = \underbrace{h_{b}(T_{b,ext}-T_{I'})}_{\phi_{b}\text{ qui intègre l'effet
de couche limite}}
  =
    \underbrace{h'_{imp,ext}(T_{imp,ext}-T_{b,ext})}_{\phi_{\text{\it réel
imposé}}}
\end{equation}

Avec pour les coefficients d'échange~:
\begin{equation}
\left\{\begin{array}{lll}
h_{int}&=&\displaystyle\frac{K}{C_v\,\overline{I'F}}\\
h_b&=&\displaystyle\frac{\phi_b}{T_{b,ext}-T_{I'}}=\frac{\rho\,C_p\,u_k}{T^+_{I'}}
\end{array}\right.
\end{equation}

On tire $T_{b,ext}$
de la seconde partie de l'égalité~(\ref{Cfbl_Cfxtcl_eq_conservation_flux_cfxtcl})
traduisant la conservation du flux~:
\begin{equation}
\displaystyle T_{b,ext} = \frac{h'_{imp,ext}\,T_{imp,ext}+h_b\,T_{I'}}{h_b+h'_{imp,ext}}
\end{equation}

En utilisant cette valeur et la première partie de l'équation de conservation
du flux~(\ref{Cfbl_Cfxtcl_eq_conservation_flux_cfxtcl}), on obtient~:
\begin{equation}
e_{b,int} = \frac{h_b\,h'_{imp,ext}}{h_{int}\,(h_b+h'_{imp,ext})}\,(T_{imp,ext}-T_{I'})+e_{I'}
\end{equation}

On utilise alors
$\displaystyle T_{I'}=\frac{1}{C_v}\left(e_{I'}-\frac{1}{2}u^2_{i}-\varepsilon_{sup,i}\right)$ pour
écrire (sans reconstruction pour la vitesse et $\varepsilon_{sup}$)~:
\begin{equation}
\displaystyle e_{b,int} =
\frac{ -\frac{h_b\,h'_{imp,ext}}{C_v}+h_{int}\,(h_b+h'_{imp,ext}) }
     { h_{int}\,(h_b+h'_{imp,ext}) } \,e_{I'}
+\frac{h_b\,h'_{imp,ext}}{h_{int}\,(h_b+h'_{imp,ext})}\,
  \left(T_{imp,ext}+\frac{\frac{1}{2}u^2_{i}+\varepsilon_{sup,i}}{C_v}\right)
\end{equation}


Et on a donc, avec $\displaystyle h'_r=\frac{h_{int}}{\frac{h_b}{C_v}}$~:
\begin{equation}
\displaystyle e_{b,int} =
\underbrace{\frac{ h'_{imp,ext} }{ C_v\,h_{int}+h'_r\,h'_{imp,ext} }\,
  \left(C_v\,T_{imp,ext}+\frac{1}{2}u^2_{i}+\varepsilon_{sup,i}\right)}_{A_b}
+\underbrace{\frac{ C_v\,h_{int}+h'_{imp,ext}(h'_r-1) }{ C_v\,h_{int}+h'_r\,h'_{imp,ext} }}_{B_b}\,e_{I'}
\end{equation}

Avec ces notations, $h_b$ est le coefficient habituellement calculé pour la
température.

Le coefficient $h'_{imp,ext}$ est le coefficient d'échange externe qui est
imposé pour la température\footnote{Le coefficient $h'_{imp,ext}$
est utile pour les cas où l'on
souhaite relaxer la condition à la limite~:
pour la température, cela correspond à imposer une valeur sur la face
externe d'une paroi unidimensionnelle idéale, sans inertie,
caractérisée par un simple coefficient d'échange.}.
Pour obtenir l'équivalent dimensionnel de $h'_{imp,ext}$ pour l'énergie,
il faut diviser sa valeur par $C_v$ (ce qui ne fait pas de différence dans
la majorité des cas, car il est habituellement pris infini).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Mise en \oe uvre}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%=================================
\subsection*{Introduction}
%=================================

Les conditions aux limites sont imposées par une suite de sous-programmes,
dans la mesure où l'on a cherché à rester cohérent avec la structure
standard de \CS.

Dans \fort{ppprcl} (appelé par \fort{precli}), on initialise les tableaux
avant le calcul des conditions aux limites~:
\begin{itemize}
\item \var{IZFPPP} (numéro de zone, inutilisé, fixé à zéro),
\item \var{IA(IIFBRU)} (repérage des faces de bord pour
lesquelles on applique un schéma de Rusanov~: initialisé à zéro,
on imposera la valeur 1 dans \fort{cfrusb} pour les faces auxquelles on applique le schéma
de Rusanov)
\item \var{IA(IIFBET)} (repérage des faces de paroi à température ou
à flux thermique imposé~: initialisé à 0, on imposera la valeur 1
dans \fort{cfxtcl} lorsque la température ou le flux est imposé),
\item \var{RCODCL(*,*,1)} (initialisé à \var{-RINFIN} en prévision
du traitement des sorties réentrantes pour lesquelles l'utilisateur
aurait fourni une valeur à imposer en Dirichlet),
\item flux convectifs de bord pour la quantité de mouvement et l'énergie
(initialisés à zéro).
\end{itemize}


\bigskip
Les types de frontière (\var{ITYPFB}) et les valeurs nécessaires
(\var{ICODCL}, \var{RCODCL}) sont imposés par l'utilisateur dans \fort{uscfcl}.

On convertit ensuite ces données dans \fort{cs\_boundary\_conditions} pour qu'elles
soient directement utilisables lors du calcul des matrices et des seconds membres.

Pour cela, \fort{cfxtcl} permet de réaliser le calcul des valeurs de bord et,
pour certaines frontières, des flux convectifs. On fait appel,
en particulier,
à \fort{uscfth} (utilisation de la thermodynamique) et à \fort{cfrusb}
(flux convectifs par le schéma de Rusanov). Lors de ces calculs, on utilise
\var{COEFA} et \var{COEFB} comme tableaux de travail (transmission de valeurs
à \fort{uscfth} en particulier) afin de renseigner \var{ICODCL} et
\var{RCODCL}.
Après \fort{cfxtcl},
le sous-programme \fort{typecl} complète quelques valeurs par défaut
de \var{ICODCL} et de \var{RCODCL}, en particulier pour les scalaires passifs.

Après \fort{cfxtcl} et \fort{typecl}, les tableaux \var{ICODCL} et \var{RCODCL}
sont complets. Ils sont utilisés dans la suite de \fort{cs\_boundary\_conditions} et en particulier
dans \fort{cs\_boundary\_conditions\_set\_coeffs\_turb} pour construire les tableaux \var{COEFA} et \var{COEFB}
(pour l'énergie, on dispose de deux couples (\var{COEFA}, \var{COEFB}) afin de
traiter les parois).

On présente ci-après les points dont l'implantation diffère
de l'approche standard. Il s'agit de
l'utilisation d'un schéma de Rusanov pour le calcul des flux convectifs
en entrée et sortie (hormis sortie supersonique)
et du mode de calcul des flux diffusifs d'énergie en paroi.
On insiste en particulier sur l'impact des conditions aux limites
sur la construction des seconds membres de l'équation de la quantité
de mouvement et de l'équation de l'énergie (\fort{cfqdmv} et \fort{cfener}).

%=================================
\subsection*{Flux de Rusanov pour le calcul des flux convectifs en entrée et sortie}
%=================================

Le schéma de Rusanov est utilisé pour calculer des flux convectifs de bord
(masse, quantité de mouvement et énergie) aux entrées et des sorties
de type IESICF, ISOPCF, IERUCF, IEQHCF.

La gestion des conditions aux limites est différente de celle adoptée
classiquement dans \CS, bien que l'on se soit efforcé de s'y conformer le
mieux possible.

En volumes finis, il faut disposer de conditions aux
limites pour trois utilisations principales au moins~:
         \begin{itemize}
        \item imposer les flux de convection,
        \item imposer les flux de diffusion,
        \item calculer les gradients pour les reconstructions.
        \end{itemize}
Dans l'approche standard de \CS, les conditions aux limites sont définies par
variable et non pas par terme discret\footnote{Par exemple, pour un scalaire
convecté et diffusé, on définit une valeur de bord unique {\it pour le scalaire}
et non pas une valeur de bord pour le {\it flux convectif} et une valeur de bord
pour le {\it flux diffusif}.}. On dispose donc, {\it pour chaque variable},
d'une valeur de bord dont devront être déduits les flux de
convection, les flux de diffusion et les gradients\footnote{Néanmoins, pour
certaines variables comme la vitesse par exemple, \CS dispose de deux valeurs
de bord (et non pas d'une seule) afin de pouvoir imposer de manière
indépendante le gradient normal et le flux de diffusion.}.
Ici, avec l'utilisation d'un schéma de
Rusanov, dans lequel le flux convectif est traité dans son ensemble,
il est impératif
de disposer d'un moyen d'imposer directement sa valeur au bord\footnote{Il
serait possible de calculer une valeur de bord fictive des variables d'état qui
permette de retrouver le flux convectif calculé par le schéma de Rusanov,
mais cette valeur ne permettrait pas d'obtenir
un flux de diffusion et un gradient satisfaisants.}.

Le flux convectif calculé par le schéma de Rusanov
sera ajouté directement au second membre
des équations de masse, de quantité de mouvement et d'énergie. Comme ce
flux contient, outre la contribution des termes convectifs ``usuels''
($\dive(\vect{Q})$, $\dive(\vect{u}\otimes\vect{Q})$ et
$\dive(\vect{Q}\,e)$), celle des termes en $\grad\,P$ (quantité de
mouvement) et $\dive(\vect{Q}\,\frac{P}{\rho})$
(énergie), il faut veiller à ne pas
ajouter une seconde fois les termes de bord issus de  $\grad\,P$ et de
$\dive(\vect{Q}\,\frac{P}{\rho})$
au second membre des équations de quantité de
mouvement et d'énergie.


Pour la masse, le flux convectif calculé par le schéma de Rusanov
définit simplement le flux de masse au bord
(\var{PROPFB(IFAC,IPPROB(IFLUMA(ISCA(IENERG))))}).

Pour la quantité de mouvement, le flux convectif calculé par le schéma de
Rusanov  est stocké dans les tableaux
\var{PROPFB(IFAC,IPPROB(IFBRHU))}, \var{PROPFB(IFAC,IPPROB(IFBRHV))} et
\var{PROPFB(IFAC,IPPROB(IFBRHW))}. Il est ensuite ajouté au second membre de
l'équation directement dans \fort{cfqdmv} (boucle sur les faces de bord).
Comme ce flux contient la contribution du terme convectif usuel
$\divv(\vect{u}\otimes\vect{Q})$, il ne faut pas l'ajouter dans
le sous-programme \fort{cfbsc2}.
De plus, le flux convectif calculé par le schéma de Rusanov
contient la contribution du
gradient de pression. Or, le gradient de pression est calculé dans
\fort{cfqdmv} au moyen de \fort{grdcel} et ajouté au second membre
sous forme de contribution volumique (par cellule)~: il faut donc retirer
la contribution des faces de bord auxquelles est appliqué le schéma de
Rusanov, pour ne pas la compter deux fois (cette opération est réalisée
dans \fort{cfqdmv}).

Pour l'énergie, le flux convectif calculé par le schéma de
Rusanov est stocké dans le tableau
\var{PROPFB(IFAC,IPPROB(IFBENE))}. Pour les faces auxquelles n'est pas
appliqué le schéma de Rusanov, on ajoute la contribution
du terme de transport de pression $\dive(\vect{Q}\,\frac{P}{\rho})$
au second membre de l'équation dans \fort{cfener}
et on complète le second membre dans \fort{cfbsc2} avec la contribution du
terme convectif usuel $\dive(\vect{Q}\,e)$. Pour les faces auxquelles est
appliqué le schéma de Rusanov, on ajoute directement le flux de Rusanov au second
membre de l'équation dans \fort{cfener}, en lieu et place de la contribution
du terme de transport de pression et l'on prend garde de ne pas
comptabiliser une seconde fois le flux convectif usuel
$\divv(\vect{Q}\,e)$ dans le sous-programme \fort{cfbsc2}.

C'est l'indicateur \var{IA(IIFBRU)}
(renseigné dans \fort{cfrusb}) qui permet, dans \fort{cfbsc2},
\fort{cfqdmv} et \fort{cfener},
de repérer les faces de bord pour lesquelles on a calculé
un flux convectif avec le schéma de Rusanov.


%=================================
\subsection*{Flux diffusif d'énergie}
%=================================

%---------------------------------
\subsubsection*{Introduction}
%---------------------------------

Une condition doit être fournie sur toutes les frontières pour le calcul du
flux diffusif d'énergie.

Il n'y a pas lieu de
s'étendre particulièrement sur le traitement de certaines frontières.
Ainsi, aux entrées et sorties, on dispose
d'une valeur de bord (issue de la résolution du problème
de Riemann)
que l'on utilise dans la formule discrète classique donnant le
flux\footnote{Les valeurs de $u^2$ et de $\varepsilon_{sup}$ ne sont pas
reconstruites pour le calcul du gradient au bord dans
$\displaystyle\dive{\left(\frac{K}{C_v}\,\grad{(\frac{1}{2}\,u^2+\varepsilon_{sup})}\right)}$}.
La situation est simple aux symétries également, où un flux nul est imposé.

Par contre, en paroi, les conditions de température ou de flux thermique
imposé doivent être traitées avec plus d'attention, en particulier
lorsqu'une couche limite turbulente est présente.

%---------------------------------
\subsubsection*{Coexistence de deux conditions de bord}
%---------------------------------

Comme indiqué dans la partie "discrétisation",
les conditions de température ou de flux conductif
imposé en paroi se traduisent,
pour le flux d'énergie, au travers du terme
$\displaystyle\dive{\left(\frac{K}{C_v}\,\grad\,e\right)}$,
en imposant une condition de flux nul sur le terme
$\displaystyle-\dive{\left(\frac{K}{C_v}\,\grad{(\frac{1}{2}\,u^2+\varepsilon_{sup})}\right)}$.
Les faces IFAC
concernées sont repérées dans \fort{cfxtcl} par l'indicateur
\var{IA(IIFBET+IFAC-1) = 1} (qui vaut 0 sinon, initialisé
dans \fort{ppprcl}).

Sur ces faces,
on calcule une valeur de bord de l'énergie, qui, introduite dans la
formule générale de flux utilisée au bord dans \CS, permettra de retouver le
flux souhaité. La valeur de bord est une simple valeur numérique sans
signification physique et ne doit être utilisée que pour calculer le flux
diffusif.

En plus de cette valeur de bord destinée à retrouver le
flux diffusif, il est nécessaire de disposer
d'une seconde valeur de bord de l'énergie afin de pouvoir en calculer le
gradient.

Ainsi, comme pour la vitesse en $k-\varepsilon$, il est nécessaire de
disposer pour l'énergie de deux couples de coefficients
(\var{COEFA},\var{COEFB}), correspondant à deux valeurs de bord distinctes,
dont l'une est utilisée pour le calcul du flux diffusif spécifiquement.

%---------------------------------
\subsubsection*{Calcul des \var{COEFA} et \var{COEFB} pour les faces de paroi
à température imposée}
%---------------------------------

Les  faces de paroi  \var{IFAC} à température imposée sont identifées par
l'utilisateur dans \fort{uscfcl} au moyen de  l'indicateur
\var{ICODCL(IFAC,ISCA(ITEMPK))=5} (noter que
ce tableau est associé à la température).

Dans \fort{cfxtcl}, on impose alors \var{ICODCL(IFAC,ISCA(IENERG))=5} et
on calcule la quantité
$C_v\,T_{imp,ext}+\frac{1}{2}u^2_{I}+\varepsilon_{sup,I}$, que l'on
stocke dans \var{RCODCL(IFAC,ISCA(IENERG),1)} (on ne reconstruit pas les
valeurs de $u^2$ et $\varepsilon_{sup}$ au bord, cf. \S\ref{Cfbl_Cfxtcl_prg_a_traiter}).

à partir de ces valeurs de \var{ICODCL} et \var{RCODCL},
on renseigne ensuite dans \fort{cs\_boundary\_conditions\_set\_coeffs\_turb}
les tableaux de conditions aux limites  permettant le calcul du flux~:
\var{COEFA(*,ICLRTP(ISCA(IENERG),ICOEFF))} et
\var{COEFB(*,ICLRTP(ISCA(IENERG),ICOEFF))} (noter
l'indicateur \var{ICOEFF} qui renvoie aux coefficients dédiés au flux
diffusif).


%---------------------------------
\subsubsection*{Calcul des \var{COEFA} et \var{COEFB} pour les faces de paroi
à flux thermique imposé}
%---------------------------------

Les  faces de paroi  \var{IFAC} à flux thermique
imposé sont identifées par
l'utilisateur dans \fort{uscfcl} au moyen de  l'indicateur
\var{ICODCL(IFAC,ISCA(ITEMPK))=3} (noter que le tableau est
associé à la température).

Dans \fort{cfxtcl}, on impose alors \var{ICODCL(IFAC,ISCA(IENERG))=3} et
on transfère la valeur du flux de  \var{RCODCL(IFAC,ISCA(ITEMPK),3)}
à \var{RCODCL(IFAC,ISCA(IENERG),3)}.

à partir de ces valeurs de \var{ICODCL} et \var{RCODCL},
on renseigne ensuite dans \fort{cs\_boundary\_conditions} les tableaux de conditions aux limites
permettant le calcul du flux,
\var{COEFA(*,ICLRTP(ISCA(IENERG),ICOEFF))} et
\var{COEFB(*,ICLRTP(ISCA(IENERG),ICOEFF))} (noter
l'indicateur \var{ICOEFF} qui renvoie aux coefficients dédiés au flux
diffusif).

%---------------------------------
\subsubsection*{Gradient de l'énergie en paroi à température ou à flux thermique imposé}
%---------------------------------

Dans les deux cas (paroi à température ou à flux thermique imposé),
on utilise les tableaux
\var{COEFA(*,ICLRTP(ISCA(II),ICOEF))},
\var{COEFB(*,ICLRTP(ISCA(II),ICOEF))} (noter le \var{ICOEF}) pour disposer d'une
condition de flux nul pour l'énergie (avec \var{II=IENERG}) et
pour la température (avec \var{II=ITEMPK})
si un calcul de gradient est requis.

Un gradient est en particulier utile pour les reconstructions
de l'énergie sur maillage non orthogonal.
Pour la température, il s'agit d'une précaution, au cas
où l'utilisateur aurait besoin d'en calculer le gradient.

%---------------------------------
\subsubsection*{Autres frontières que les parois à température ou à flux thermique imposé}
%---------------------------------

Pour les frontières qui ne sont pas des parois à température ou
à flux thermique imposé, les conditions aux limites de l'énergie et
de la température sont complétées classiquement dans \fort{cs\_boundary\_conditions} selon
les choix faits dans \fort{cfxtcl} pour \var{ICODCL} et \var{RCODCL}.

En particulier,
dans le cas de conditions de Dirichlet sur l'énergie (entrées, sorties), les
deux jeux de conditions aux limites sont identiques (tableaux
\var{COEFA}, \var{COEFB} avec \var{ICOEFF} et \var{ICOEF}).

Si un flux est imposé pour l'énergie totale (condition assez rare,
l'utilisateur ne raisonnant pas,
d'ordinaire, en énergie totale), on le stocke au moyen de
\var{COEFA(*,ICLRTP(ISCA(IENERG),ICOEFF))} et
\var{COEFB(*,ICLRTP(ISCA(IENERG),ICOEFF))} (tableaux associés au flux
diffusif). Pour le gradient, une condition de flux nul est stockée
dans
\var{COEFA(*,ICLRTP(ISCA(IENERG),ICOEF))} et
\var{COEFB(*,ICLRTP(ISCA(IENERG),ICOEF))}. On peut remarquer que les deux
jeux de conditions aux limites sont identiques pour les faces de symétrie.

%---------------------------------
\subsubsection*{Impact dans \fort{cfener}}
%---------------------------------

Lors de la construction des seconds membres, dans \fort{cfener}, on utilise les
conditions aux limites stockées dans les tableaux associés au flux
diffusif
\var{COEFA(*,ICLRTP(ISCA(IENERG),ICOEFF))} et
\var{COEFB(*,ICLRTP(ISCA(IENERG),ICOEFF))} pour le terme de flux diffusif
$\displaystyle\dive{\left(\frac{K}{C_v}\,\grad\,e\right)}$
en prenant soin d'annuler la contribution de bord du terme
$\displaystyle-\dive{\left(\frac{K}{C_v}\,\grad{(\frac{1}{2}\,u^2+\varepsilon_{sup})}\right)}$
sur les faces pour lesquelles cette condition
prend les deux termes en compte, c'est-à-dire sur les faces pour lesquelles
\var{IA(IIFBET+IFAC-1) = 1}.

Pour tous les autres termes qui requièrent une valeur de bord, on utilise les
conditions aux limites que l'on a stockées au moyen des deux tableaux
\var{COEFA(*,ICLRTP(ISCA(IENERG),ICOEF))} et
\var{COEFB(*,ICLRTP(ISCA(IENERG),ICOEF))}. Ces conditions sont
donc en particulier utilisées pour le calcul du gradient de l'énergie,
lors des reconstructions sur maillage non orthogonal.


\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Points à traiter}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\label{Cfbl_Cfxtcl_prg_a_traiter}%
% propose en patch 1.2.1
%Corriger \fort{ppprcl} pour que l'indicateur
%\var{IA(IIFBET+IFAC-1)} soit
%initialisé à 0, positionné à 1 aux faces de paroi à température
%ou flux thermique imposé. Dans \fort{cfener}, lorsque l'indicateur vaudra 1,
%on ne prendra pas en compte le flux correspondant à
%$\displaystyle-\dive{\left(\frac{K}{C_v}\,\grad{(\frac{1}{2}\,u^2+\varepsilon_{sup})}\right)}$.

% propose en patch 1.2.1
%Pour l'énergie, on utilise comme diffusivité turbulente la valeur
%$\displaystyle \frac{K}{C_v}=\frac{\lambda}{C_v}+\frac{\mu_t}{\sigma_t}$.
%Par cohérence avec une équation
%portant sur la température, il serait plus logique d'utiliser
%$\displaystyle \frac{K}{C_v}=\frac{\lambda}{C_v}+\frac{C_p}{C_v}\,\frac{\mu_t}{\sigma_t}$.
%On peut temporairement utiliser le nombre de Prandtl turbulent pour prendre en compte
%le rapport $\displaystyle\frac{C_p}{C_v}$, mais il serait
%souhaitable de corriger en ce sens le calcul de \var{W1} pour \fort{cs\_face\_viscosity} dans
%le sous-programme \fort{cfener} et le calcul similaire de \var{HINT} dans
%\fort{cs\_boundary\_conditions} et \fort{cs\_boundary\_conditions\_set\_coeffs\_turb} (RAS pour les conversions en couplage avec \syrthes).

Apporter un complément de test sur une cavité fermée
sans vitesse et sans gravité, avec flux de bord ou température de bord imposée.
Il semble que le transfert d'énergie {\it via} les termes de pression génère de
fortes vitesses non physiques dans la première maille de paroi et que la
conduction thermique ne parvienne pas à établir le profil de température
recherché. Il est également possible que la condition de bord sur la pression
génère une perturbation (une extrapolation pourrait se révéler
indispensable).

Il pourrait être utile de généraliser à l'incompressible l'approche
utilisée en compressible pour unifier simplement le traitement
des sorties de type 9 et 10.

Il pourrait être utile d'étudier plus en détail l'influence de la non
orthogonalité des mailles en sortie supersonique (pas de reconstruction,
ce qui n'est pas consistant pour les flux de diffusion).

De même, il serait utile d'étudier l'influence de l'absence de
reconstruction pour la vitesse et $\varepsilon_{sup}$ dans la relation
$\displaystyle T_{I'}=\frac{1}{C_v}\left(e_{I'}-\frac{1}{2}u^2_{i}-\varepsilon_{sup,i}\right)$
utilisée pour les parois à température imposée.

Apporter un complément de documentation pour le couplage avec \syrthes (conversion
énergie température). Ce n'est pas une priorité.

Pour les thermodynamiques à $\gamma$ variable, il sera nécessaire de
modifier non
seulement \fort{uscfth} mais également \fort{cfrusb} qui doit disposer de
$\gamma$ en argument.

Pour les thermodynamiques à $C_v$ variable, il sera nécessaire de
prendre en compte un terme en $\grad\,C_v$, issu des flux diffusifs,
au second membre de l'équation de
l'énergie (on pourra cependant remarquer qu'actuellement, en incompressible,
on néglige le terme en $\grad\,C_p$ dans l'équation de l'enthalpie).
